In [97]:
import matplotlib
#matplotlib.use('Agg')
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
In [98]:
import numpy as np
#from chainconsumer import ChainConsumer
from corner import corner
In [99]:
! ls -lt /scratch/users/swmclau2/PearceMCMC/*.npy
In [100]:
from glob import glob
fnames = glob('/scratch/users/swmclau2/PearceMCMC/500_walkers_20000_steps_xigg_emu2_*.npy')
In [101]:
n_walkers = 500
n_burn = 15000
In [102]:
from itertools import islice
In [103]:
def read_chain(fname, n_walkers, n_burn)
chain_list = []
with open(fname, 'r') as f:
step_counter = 0
header = f.readline() #header
chain_pnames = header[1:].split()
while True:
next_lines = islice(f, n_walkers)
step_counter+=1
#if step_counter % 1000 == 0:
# print step_counter, len(chain_list)
if step_counter <= n_burn:
for line in next_lines:
pass
continue
next_lines = np.array([ np.fromstring(line, sep = ' ') for line in next_lines], dtype=float)
if len(next_lines) < 1:
break
chain_list.append(next_lines)
return np.vstack(chain_list)
In [ ]:
chains = [read_chain(fname, n_walkers, n_burn) for fname in fnames]
In [106]:
n_params = chain.shape[1] if len(chain.shape) > 1 else 1
In [108]:
param_names = [r'$N_{eff}$', r'$\log(M_0)$',r'$\log(M_1)$', r'$H_0$',r'$w_0$', r'$\ln(10A_s)$', r'$\Omega_c h^2$',
r'$\sigma_{\log M }$', r'$\alpha$', r'$n_s$', r'$\Omega_b h^2$']
In [109]:
hod_idxs = np.array([1, 2, 7, 8])
cosmo_idxs = np.array([0, 3, 4, 5, 6, 9, 10])
In [110]:
#hod_chain = chain[:, hod_idxs]
#cosmo_chain = chain[:, cosmo_idxs]
cosmo_chain = chain
In [111]:
chain_pnames
Out[111]:
In [112]:
hod_param_names = [r'$\log(M_0)$',r'$\log(M_1)$', r'$\sigma_{\log M }$' ,r'$\alpha$' ]
cosmo_param_names = [r'$N_{eff}$', r'$H_0$', r'$w_0$', r'$\Omega_c h^2$', r'$\ln(10A_s)$' ,r'$n_s$',r'$\Omega_b h^2$' ]
In [113]:
from pearce.mocks import cat_dict
cosmo_params = {'simname':'testbox', 'boxno': 0, 'realization':0, 'scale_factors':[1.0], 'system': 'sherlock'}
cat = cat_dict[cosmo_params['simname']](**cosmo_params)#construct the specified catalog!
In [114]:
cpv = cat._get_cosmo_param_names_vals()
cat_val_dict = {key: val for key, val in zip(cpv[0], cpv[1])}
In [115]:
#cosmo_true_vals = [3.7,70.7317,-1.13151,0.12283, 3.11395, 0.953515, 0.021762]
cosmo_true_vals = [cat_val_dict[pn] for pn in chain_pnames if pn in cat_val_dict]
In [116]:
print cosmo_true_vals
In [117]:
emulation_point = [('logM0', 14.0), ('sigma_logM', 0.2),
('alpha', 1.083),('logM1', 13.7)]
In [ ]:
hod_true_vals = [14.0, 13.7, 0.2, 1.083]
In [ ]:
corner(cosmo_chain, labels=cosmo_param_names,
quantiles=[0.13, 0.5, 0.86],
truths = cosmo_true_vals,
show_titles=True, title_kwargs={"fontsize": 12},
plot_datapoints = True, plot_density = True);
In [ ]: